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ABSTRACT 



Efforts to detect gravitational waves by timing an array of pulsars have focused tra- 
ditionally on stationary gravitational waves: e.g., stochastic or periodic signals. Grav- 
itational wave bursts — signals whose duration is much shorter than the observation 
period — will also arise in the pulsar timing array waveband. Sources that give rise 
to detectable bursts include the formation or coalescence of supermassive black holes 
(SMBHs), the periapsis passage of compact objects in highly elliptic or unbound orbits 
about a SMBH, or cusps on cosmic strings. Here we describe how pulsar timing array 
data may be analyzed to detect and characterize these bursts. Our analysis addresses, 
in a mutually consistent manner, a hierarchy of three questions: i) What are the odds 
that a dataset includes the signal from a gravitational wave burst? u) Assuming the 
presence of a burst, what is the direction to its source? and Hi) Assuming the burst 
propagation direction, what is the burst waveform's time dependence in each of its po- 
larization states? Applying our analysis to synthetic data sets we find that we can detect 
gravitational waves even when the radiation is too weak to either localize the source of 
infer the waveform, and detect and localize sources even when the radiation amplitude is 
too weak to permit the waveform to be determined. While the context of our discussion 
is gravitational wave detection via pulsar timing arrays, the analysis itself is directly 
applicable to gravitational wave detection using either ground or space-based detector 
data. 

Subject headings: methods: data analysis — methods: statistical — gravitational waves 



Introduction 



It has been just over thirty years since Sazhin (1978) and Detweiler (1979) showed how gravi- 



tational waves could be detected by correlating the timing residuals of a collection of pulsars, and 
twenty years since Foster & Backer (1990) proposed using a collection of pulsars — i.e., a pulsar 
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timing array — to achieve greater sensitivity. Over the ensuing years telescope collecting area has 
increased, antenna temperature has decreased, pulsar timing electronics and methodology has im- 
proved, and pulsars with exceptionally low intrinsic timing noise have been discovered. As a result 
of these advances the near-future detection of a stochastic gravitational wave signal through pulsar 
timing observations is a strong possibility. Analyses aimed at detecting gravitational waves using 
pulsar timing array observations have traditionally focused on stationary signals (i.e., stochastic or 
periodic gravitational waves). More generally, analyses aimed at detecting gravitational waves have 
merged the questions of detection and characterization, overlooking the possibility of detecting a 
signal that is too weak to be characterized. Here we describe how pulsar timing array data may be 
analyzed to search for gravitational wave bursts, demonstrating that i) pulsar timing array data 
is sufficiently rich to allow the detection of gravitational wave bursts, the localization of the burst 
source, and the time-dependent waveform of the radiation in its (two) polarization states; and, ii) 
that gravitational wave signals too weak to be characterized, or too weak to allow their source to 
be localized, may still be strong enough to be unambiguously detected. 

The first detections of gravitational waves will be important for confirming their existence 
and testing whether general relativity correctly predicts their properties (e.g., polarization modes, 
propagation speed). Of perhaps greater long-term significance will be the use of gravitational 
waves as a tool of observational astronomy that gives us direct insight into phenomena that we can 



now observe only indirectly, if at all. For example, Jaffe & Backer (2003); Wyithe & Loeb (2003) 



and Jenet et al. (2006) have shown that root-mean-square (rms) amplitude of a stochastic signal 



arising from the confusion limit of a large number of supermassive black hole binary coalescences is 
within an order of magnitude of the current sensitivity of the most advanced pulsar timing array. 
Since the signals that contribute to this background arise from a population of discrete sources 
distributed throughout space we quite reasonably expect that some of the individual sources may be 



observable as gravitational wave bursts rising above this background. Indeed, recent work by Sesana 
et al. (2008) shows that at frequencies greater than a few times 10 -8 nanohertz the gravitational 



"background" arising from supermassive binary black hole coalescense should be dominated by a 
few bright sources. Other potential burst gravitational wave sources in the pulsar timing array 
band include cosmic (super) string cusps and kinks ( Damour Sz Vilenkin|2001~ Siemens et al.| [2007 



Leblond et al. 2009) 



Pulsar timing array observations are sensitive to gravitational waves of periods ranging from 
the interval between timing observations (days to months) and the duration of the observational 
data sets (years). The corresponding wavelengths are much greater than those explored in existing 
or proposed human-built ground or space-based detectors. Ground-based detectors, whether of 
the acoustic (Astone et al. 2010) or interferometric variety ( Accadia et al.pOlO Riles et al. 2010), 
are currently sensitive to waves in the ~ 100 Hz - 1 kHz band with proposed advances opening- 



up the 10-100 Hz band (Smith & the LIGO Scientific Collaboration 2009; Kuroda & the LCGT 



Collaboration|2006 ) . Space-based detectors, which have been the subject of extensive design studies 
over the last thirty years, would be sensitive to gravitational waves in the ~ 3 x 10" 5 Hz - 10 Hz 
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band ( |Stebbins|200^|Jennrich|2009|pCawamura et al.|2008[ ). Over this broad band — 10" 9 - 10 3 Hz 
- the scale and character of the sources varies dramatically: e.g., ground-based detectors will be 
sensitive to gravitational waves from neutron star or solar mass black hole binaries, supernovae 
and gamma-ray burst progenitors; space-based detectors to gravitational waves from white dwarf 
binaries, stellar disruptions about intermediate mass black holes and the inspiral of solar mass 

compact objects or intermediate mass black holes about 10 4 ' 5 10 7,5 M black holes; and pulsar 

timing arrays to the formation, interaction and evolution of supermassive black holes. Pulsar timing 
array observations thus offer their own, unique perspective on the gravitational wave universe, 
distinct from that provided by the either ground- and space-based detectors. 

Analyses aimed at detecting gravitational waves using pulsar timing array observations have 
traditionally focused on stationary signals: i.e., an isotropic stochastic gravitational wave back- 



ground (jHellings & Downs 1983 


McHugh et al. 1996 Thorsett k Dewey 


1996 


Lommen||2001 


Lommen et al. 


2003 Jenet et al. 


2005 Jenet et al. 2006 Demorest |2007 Hobbs et al.||2008[ van 


Haasteren et al 


. 2009 Anholm et al. 2009 


) or gravitational waves from discrete periodic sources 


(Lommen & Backer|2001 Jenet et al.|;2004 


Jenet et al.|2005a|b). More recent work 


van Haasteren 


k Levin|2009 ) has investigated the detection of gravitational wave "memory" 


(Christodoulou|1991) 



associated with sources that radiate a significant amount of energy in gravitational waves (Wiseman 



k Will|[l99T| ) or that become unbound ( |Thorne|[T992| ) . 



Gravitational wave detection using a pulsar timing array begins with the collection of timing 
residuals from the several array pulsars. These timing residuals are the difference between the 
expected pulse arrival times (taking into account all non-gravitational- wave propagation effects) 
and the actual pulse arrival times at each pulsar observational epoch. For pulsars used in current 
timing arrays the timing precision is in the 50 ns - 5 /is range. In ^2]we summarize how these timing 
residuals reflect the passage of a plane gravitational wave through the pulsar-Earth baseline. In £j3] 
we describe our analysis for gravitational wave bursts, which takes advantage of the correlation of 
the timing residuals measured for different pulsars. In ^4] we demonstrate the effectiveness of the 
analysis by applying it to simulated data arising from a thirty pulsar timing array and including a 
gravitational wave burst such as would be expected from a parabolic encounter of two supermassive 
black holes. Finally, in §[5] we summarize our findings and describe planned future work. 



2. Pulsar timing response to the passage of a gravitational wave burst 

2.1. Introduction 

A pulsar timing array dataset consists of a collection of pulsar "time of arrival" , or TOA, mea- 
surements for pulses of the individual pulsars that comprise the array. The arrival time observations 
are made for each pulsar over a period years, with successive pulse arrival time observations for each 
array pulsar made anywhere from days to months apart. The TOA measurements are compared 
to predicted arrival times based on timing models for the individual pulsars, which includes all 
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non-gravitational-wave effects that affect the arrival times. The difference between the observed 
and expected pulse arrival times are referred to as timing residuals, which are then presumed to 
consist of timing noise and gravitational wave effects. Evidence for gravitational waves is sought 
in the timing residuals [j] In this section we calculate the contribution to pulse arrival times owing 
to a passing plane gravitational wave burst. 



2.2. Gravitational waves 

Denote the perturbative plane gravitational wave, expressed in transverse-traceless (TT) gauge 
( |Misner et aT|l973| ), as 



h(t,x) = h + {t-k-x)e {+) (k) + h x {t-k-x)e^ x) (k) (2-1) 

where k is the plane wave propagation direction and e^ + ) and e^ x ) are the two independent gravi- 
tational wave polarization basis tensors, 

Jm ( + ) _ Im (X) _ cy /r, r, \ 

e (+) e /m - e (x) e /m - Z 

e l ™)k m = e l j™jk m = efye ln ? = 0. (2-2b) 

Locating the coordinate system origin at the solar system barycentre consider a pulsar at spatial 
rest located at x p , 

x p {t) = Ln (2-3) 
where n is the unit vector in the direction of the pulsar and L is pulsar's distance. 



2.3. Timing residuals 

Focus attention on the electromagnetic field associated with the pulsed emission of a pulsar 
and denote the fields phase, at the pulsar, as (f>o(t). We are interested in the time-dependent phase 
(j)(t) of the electromagnetic field associated with the pulsed emission measured at an Earth-based 
radio telescope, which we write as 

<j)(t) = Mt ~ L - r {t) - TGW (t)} (2-4a) 



1 The procedure of fitting the timing model to the pulsar arrival time measurements for gravitational wave analysis 
has the unfortunate side-effect of "fitting out" any gravitational wave contributions that have the form of other timing 
model effects. We address this point directly in the conclusions. 
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where 



TO 



TGW 



^ Corrections owing exclusively to the spatial motion of the Earth 

within the solar system, the solar system with respect to the pulsar, 
i and electromagnetic wave propagation the interstellar medium 

Corrections owing exclusively to h(t, x) 



(2-4b) 
(2-4c) 



(Note that we work in units where c = G = 1.) In the absence of gravitational waves tgw vanishes 
and the front 4>o{t) arrives at Earth at time £©(i) = t + L + tq{£) . In the presence of a gravitational 
wave signal the phase front arrives at time te(*) + t gw(*); thus, tqw is the gravitational wave 



timing residual. Following Finn (2009) Eqs. (3.26) and (3.12e), the arrival time correction tqw(^) 
is 



tgw (t) = ■ 

where 

It is convenient to introduce f A (u), 



1 



n l n m 



Ira H (x) 



h A (t - (1 + k jn j ){L - A)) dX 



dfA 
du 



h A (u), 



and rewrite equation 2-5b using Ja(u) as follows: 

fA(t) 



f A (t - (1 + kjnP)L) 



1 + kjni 



1 + kj-nj 



(2-5a) 



(2-5b) 



(2-6) 



(2-7) 



The contribution proportional to f A (t) is colloquially referred to as the "Earth" term; similarly, the 
contribution proportional f A \ t — (1 + A: m n m )L^ is referred to as the "Pulsar" term. The Pulsar 
term is of central importance when pulsar timing data is used to bound the strength of a stochastic 
gravitational wave background ( Jenet et al.|2005 ); however, as we show below, only the Earth term 
is important when our goal is to use pulsar timing data to detect gravitational wave bursts. 



2.4. Discussion 



At this point it is worth noting several properties of the timing residual tqw- 



2.4-1 ■ Burst detection involves only the Earth term 



As shown in equation (2-7) the gravitational wave induced timing residuals for any pulsar may 
be written as the difference of two functions, each of which is an integral of /i+, x (t, x). These two 
functions are identical, except that one is displaced in time by an amount L{\ + k m n m ) with respect 
to the other. Correspondingly, 
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When timing residual measurements from an array of pulsars are available the first evidence 
for the passage of a gravitational wave burst will appear simultaneously in all observed resid- 
uals; and 

As long as the burst duration AT and the observation duration T are less than (l + k m h m )L/c 
only the Earth term contributes to the correlated timing residuals in the pulsar timing arrayj^] 



When searching for gravitational wave bursts we can thus ignore the pulsar term except for sources 
within an angle 



9 P < cos 



-i 



1 



at" 



2.5 deg 



"ATlkpc 
lyr L 



.1/2 



(2- 



of pulsar p. 



2.4-2. Timing residuals in a pulsar network are sensitive to gravitational wave polarization 
The timing residual tqw is a linear combination of the integrals of the two polarizations of the 



waveform H(a) : i- e -> we may rewrite 



2-5 



as 



r G w(t) = ~\ [F + n {+) + T x ^ (x) ] , (2-9) 

where 

= n l h m e\^{k). (2-10) 

The timing residual correlations of PTA pulsars take a form that depends on the pulsar locations 
and the gravitational wave polarization. When the wave propagation direction k is known the 
measured timing residuals of two appropriately chosen pulsars is sufficient to measure separately 
the radiation in each of the two gravitational wave polarization states. 



2.4-3. Timing residuals in a pulsar timing array are sensitive to wave propagation direction. 

(A) 

The polarization tensors e)' are orthogonal to the wave propagation direction k; correspond- 
ingly, the relative contribution of the H(a) to the timing residual tqw for a given pulsar depends 
on the gravitational wave propagation direction through the F^ A \ In addition, the overall ampli- 
tude of the timing residual for any particular pulsar depends on the wave propagation direction 
through the additional factor (1 + fc m n™) _1 . Observation of the timing residuals in three pulsars, 



2 Other bursts, having interacted with individual pulsars at much earlier times (thousands of years) will contribute 
to the timing noise of individual pulsars. These contributions will not be correlated among the pulsars in the timing 
array over the human observational timescale (decades). 
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with appropriately chosen lines-of-sight from Earth, are thus sufficient to measure the radiation 
propagation direction. 



Combining the insights of subsections 2.4.1 2.4.2 and 2.4.3 we see that a pulsar timing array 
of five or more pulsars has, in principle, sufficient information to fully characterize a passing grav- 
itational wave burst. In the following section we describe the statistical methodology by which we 
can infer k and h +)X (t) at, e.g., the solar system barycentre from the measured timing residuals in 
a pulsar timing array. 



2.4-4- Pulsar timing residuals are larger for longer bursts than for shorter bursts 
The gravitational wave induced timing residual associated with any particular pulsar is pro- 



portional to the integral of hij(t) over time (see eq. 2-5b). This leads to an important point: for 
fixed strain amplitude and waveform "shape", the timing residuals associated with bursts have 
magnitudes proportional to the burst duration. This is very different than is the case with ground- 
based gravitational wave detectors (e.g., the Laser Interferometer Gravitational Wave Observatory 
(LIGO) QSaulson|[i"994"l )) or the proposed space-based detector LISA, where the measured quantity 



responds directly to the gravitational wave strain. The difference arises because the gravitational 
wave signal band of interest for ground- and space-based detectors has wavelengths greater than 
the detector size, while the band of interest for pulsar timing array measurements has wavelengths 
much smaller than the detector size (i.e., the pulsar-Earth baseline distance)]^] 



3. Statistical Methodology 

3.1. Framing the questions 

Our goal is three- fold. First, ascertain the odds that the particular data set d includes a 
contribution characteristic of a passing gravitational wave burst; second, assuming that is so, de- 
termine the probability that the contribution is characteristic of a wave propagating in the direction 
k, and, finally, assuming the contribution is characteristic of a burst propagating in direction k, 
determine the probability that the contribution is characteristic of a waveform at Earth described 
by h = h + {t — k ■ x)e+(k) + h x (t — k ■ x)e x (k) for functions h + and h x ■ 

While actual analysis might address these questions in the order given above it is pedagogically 
simpler and more instructive to approach them in the opposite order, which we do in the three 
subsections that follow. 



3 For LISA the detector bandwidth does extend to wave frequencies a few times greater than the round-trip travel 
along the 5 x 10 6 km arm baseline. This effect of greater sensitivity at longer periods is apparent in the high-frequency 
part of LISA's response function. 
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3.2. Inferring h 

Given timing residual observations d from an array of pulsars that include a contribution from 
a plane gravitational wave propagating past Earth in direction k, what is the probability ph that 
the wave is described by h? 

The desired probability density depends on the response of the pulsar network to incident 
gravitational waves, the statistical properties of the measurement and intrinsic pulsar timing noise, 
and the assumed direction of wave propagation: 

f / probability that gravitational wave burst is described by the wave h \ » . 

Ph(h\k,l) = - . . (3-la) 

\ propagating in direction k, and other, unenumerated assumptions X J 

Exploiting the Bayes' Theorem, the probability density ph can be expressed in terms of the nor- 
malized likelihood A, an a priori probability density q\ t that expresses expectations regarding h, 
and a normalization constant Z^, often referred to as the evidence: 




A(d\h,k,X)g h (h\k,l) 



Z h (d\k,l) 

Probability of observing TOA residuals d given \ 
gravitational wave h propagating in direction k J 

a priori probability density that h describes the 
gravitational wave burst propagating in direction k 



(3-2a) 
(3-2b) 
(3-2c) 



(Th+dTK A(d|h,fc,I)g fc (h|fc,X) 

p d (d\k,l) (3-2d) 

(Probability of observing d assuming the presence of \ 
gravitational wave burst h propagating in direction k J 



(In equation 3-2d the integral is over all possible values of the waveform /i + and h x at the n sample 



times.) We discuss each of these terms in more detail below. 



3.2.1. The Likelihood A 

Focus attention on pulsar j, whose measured timing residuals are represented as the time- 
series vector dj. These residuals are the sum of of measurement noise, intrinsic pulsar timing noise, 
scintillation and other propagation noises nj and the pulse arrival time disturbance owing to the 
passing gravitational wave. Denoting by Rj the the timing residual response function for pulsar j 
the net timing residual measured for pulsar j is 



dj = iij + Rjh. 



(3-3) 
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The noise associated with individual pulsar timing residual observations are generally well-modeled 
as Gaussian distributed with zero mean; correspondingly, the noise associated with collection of 
observations dj are described by a zero-mean multi-variate Gaussian. Denoting the noise auto- 
correlation for pulsar j as Cj(t[ — t m ) write the probability density of observing residuals dj in the 
timing data of pulsar j as 

Aj(dj\h, k,X) = N(dj - Rjh|Cj-), (3-4a) 
where Cj is the noise auto-corrleation in detector j and 

/ (multivariate) normal distribution for zero \ /„ n 

^(x|C) = V ' . . . (3-4b) 

\ mean random deviate x given co-variance L / 



exp [-|x T C- 1 x] 



(27r) dimx det||C| 



(3-4c) 



Recall that the noise covariance C has elements 

C jk =< n(tj)n(t k ) > (3-5) 

where n(t) is the noise at time t and <> denotes an ensemble average over the noise. Expressed 
as a function of r = t k — tj, C(r) is the noise auto-correlation function, which is just the cosine- 
transform of, and thus entirely equivalent to, the noise power spectral density ( ]Kittel|1958| ) . White, 
pink, red, or more complex noise timing noise spectra are thus equally well described by Equation 

O- 

Now assume that the timing noise associated with the observations of the n p different 
pulsars are uncorrelated. Under this assumption the probability density of a set of timing residuals 
d, consisting of residuals dj from each pulsar j in the network, is 

A(d|h,jfe,X) =Y[A j (d j \h,k,l) (3-6a) 
j'=i 

=iV(d-Rh|C). (3-6b) 



3.2.2. The. prior q h 

The a priori probability density describes our expectations, before interpreting the obser- 
vations d, regarding the gravitational wave burst h. It is often the case that discussions of priors 
like these are more heated and intense than is warranted by the difference any reasonable choice 
makes to the final result. To understand how this is so it is worthwhile to return for a moment 



to Equation (3-2). The probability is the product of two /i-dependent terms, A and q^. All of 
the data dependence is encapsulated in the likelihood A; i.e., the prior is independent of the 
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observations d. When the data are conclusive A is more sharply peaked than qh and the dependence 
of ph on h is dominated by the data-dependent term A. In this case the prior qh is approximately 
constant over the volume of h where ph is large and the particular choice of prior is unimportant. 
On the other hand, when the data are inconclusive the dependence of ph on h is dominated by 
the prior and the structure of A is unimportant. As long as the prior, viewed by itself, does not 
reflect an overly strong set of expectations about h it will not matter what particular form it takes 
except at the margins where the observations are suggestive but not conclusive. With this in mind 
we consider the basic assumptions we make regarding a gravitational wave burst and how those are 
represented in qh- 

To begin, we make no assumption that the nature of the burst should be correlated with its 
direction of propagation; i.e., we drop the dependence of qh on k: 

q h (h\k,l) = q h (h\l). (3-7) 

We also assume that there is no a priori correlation between the two dynamically independent 
polarization states, in which case 

q h (h\l) = q + (h + \l)q x (h x \T) (3-8) 

where the + and x subscripts denote any two orthogonal polarization states. Since the resolution of 
a gravitational wave into orthogonal polarization states is determined only up to a rotation about 
the propagation direction it must be the case that g+ and q x are the same function qo of their 
arguments: i.e., 

q h (h\l) = q + (h+\l)q x (h x \X) = q (h + \l)q (h x \l) (3-9) 

Now suppose we represent the gravitational waveform h by the values of h + and h x at the solar 
system barycentre sampled at nh times ty. 

h +d = h+(tj) (3-10a) 
h x ,j = h x {tj). (3-10b) 

Assuming that the product /i + (i)/i + (i + r) (similarly h x {t)h x {t + r)) vanishes for r / when 



averaged over the ensemble of all possible waveforms h + (h x ), Summerscales et al. (2008) showed 
that we obtain a functional equation for go whose solution is 



q (h\a,l) = N(h\al) (3-lla) 

dimh 7 2 \ 

< 3 -" b > 



/„ 2\dim/i 
27TCT ) 



-1/2 / l^h 2 ' 



cxp 

where a is an undetermined constant and I denotes the appropriately dimensioned identity matrix. 



2 ^ a 2 

k=l 



Our minimal assumptions thus fix the prior q^ up to two undetermined constants cr+ and a x : 

q h (h\a+,a x ,l) = N(h + \a+I)N(h x \a x l). (3-12) 
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In the statistics literature the new constants <7+ and cr x are referred to as hyperparameters ( Gelman 



et al. 2004, Chapter 5). Often times the hyperparameters may have a physical interpretation 
that allows their values to be set, or a priori probability distributions (hyperpriors) selected to 
describe them, in which case the hyperparameters are treated on par with the other problem 
parameters. In our case there is a natural interpretation of cr+ jX as the root-mean-square amplitude 
of the gravitational wave burst. This interpretation is not sufficient to determine c +iX a priori or 
determine an a priori probability density over the <7+ jX . This situation is not at all uncommon. 
Several methods have been suggested and investigated for the treatment of hyperparameters in this 



case 


Galatsanos & Katsaggelos 


1992 


Thompson & Kay 


1993 


Keren & Werman 


1996 


MacKay 


1996 


Galatsanos et al. 


1998 


MacKay 


1999; 


Cawley & Talbot 


2007 


). Comparative studies suggest 



that the best treatment assigns to the hyperparameters those values that optimize the evidence Z^ 
regarded as a function of the hyperparameters ( MacKay]|1996 1999; Molina et al.||1999 ). We adopt 
this procedure here. 



The evidence Zh is the integral of Aq^ over all h x (see Eq. 3-2d). Since all of the probability 
densities that arise in our problem are normal distributions the evidence may be computed in closed 
form. Combining Eq. |3-2d| with Eqs. 3-4, 3-6 and 3-11 and completing the square in the exponential 
we obtain: 



Z h (d\a + ,a x ,l) 



exp [- 



i d l exp 



' x d\ T A~ x 



R T C x d 



v / (2vr) dimd det||C| 



A 1 1 2 dim h + 2 dim h > 
\\cr, CT V 



(3-13a) 



det 



where 



A=( a + l+ ° 2 Ur^R 
V a- 2 l x ) 



(3-13b) 



and I +iX represent the appropriately dimensioned unity matrices. 



3.2.3. The probability density ph 
To summarize, the posterior probability density pt is given by 

p h (h\k,d, a + ,a x ,T) - 

where ho satisfies 

Ah n = R T C _1 d 



det || A| | 

\dim h 



(2vr) c 



^(h-h ) T A(h-h ) 



(3-14a) 



(3-14b) 



with A given by Equation 3- 13b 



The reader may note that Equation 3-14b for ho bears a superficial resemblance to a "(reg- 
ularized) least squares" estimate for the incident wave. This resemblance is an accident of the 
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notation. The operator A that appears in Equation 3-14b would be a constant in a least squares, 
or regularized least squares, analysis. Here, however, the regularization constants a + and a x that 
appear in A get their values through the optimization of the evidence Z^, which involves both A 
and the observations d. Equation 3-14b for ho must be solved simultaneously with the optimization 
of Zfa leading to cr_|_, cr x and ho that differ from any "least squares" analysis. Finally, the principal 
result of our analysis — i.e., Equation 3-14a for p^ — would never arise from a least squares (or a 
maximum likelihood) analysis. 

As is apparent from Eq. |3-14 ho is the waveform that maximizes the probability density ph- 
As such it is naturally the "best guess" for h. The availability of the overall probability density ph 
gives us the opportunity to say and do much more. With p^ comes the ability to characterize the 
certainty should assign to this inferrence and, in general, the ability to propagate errors through 
any inferences that depend on our estimate for h. (See, for example, Bondarescu et al. (2010), 
where pn is used to estimate the uncertainty in the gravitational wave Stokes Parameters.) 



3.3. Inferring the wave propagation direction k 

Given timing residual observations d from an array of pulsars that is assumed to include the 
signal from a plane gravitational wave propagating past Earth in an unknown direction, what is 
the probability that the wave is propagating in direction k? 

The desired probability depends on the response of the pulsar network to incident waves and 
the statistical properties of the measurement and intrinsic timing noise: 

probability that burst is propagating 
Pk(k\d,X) = [ in the direction k, given data d and | (3-15) 
other, unenumerated assumptions I 

Exploiting Bayes' Theorem the probability density pk can be expressed in terms of pd, an a priori 
probability density that expresses our assumptions regarding k, and a new normalization constant, 
referred to as the evidence for k: 

p k (k\d,X) = Z^\d\l)p d (d\k,I)qk(k\I) (3-16a) 

where 

/a priori probability that the gravitational \ . 

q k {k\T) = . t \ (3-16b) 
\ wave burst is propagating in direction k J 

Z-\d\Z) = [ d 2 n kPd (d\k,l)q k (k\l) (3-16c) 
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Noting that 

P d (d\k,X) = J d n h + d n h x Pdh {d,h\k,X) (3-16d) 

= J d n h + d n h x Ph (h\d,k,l)Z h (d\k,l) (3-16e) 

= Z h (d\k,X) (3-16f) 



we find 



where 



Pk(k\d,l) = ^^ V^Z) (3-16g) 

Z k {d\X) = f d 2 fl k Z h (d\k,X)q k (k\X) (3-16h) 

A non-controversial choice of prior q k arises from assuming that we have no a priori reason to 
believe that gravitational wave bursts are propagating in any direction preferentially, in which case 
q k is uniform on the sphere (i.e., qk(k) = (47r) _1 ). In that case q(k\X) is independent of h and we 
have 

tl\A t\ 1 z h(d\k,X) 

Pkm > I) = ^^(dW ( } 

Z k (d\X) = d 2 n k Z h (d\k,l) (3-17b) 



with Zh given by equation (3-2d). 



3.4. Inferring the odds that a gravitational wave is present 

3.4-1- Model comparison and the Bayes Factor 

Given timing residual observations d from an array of pulsars, what odds should we give that 
a plane gravitational wave was incident on Earth over the period of the observation? 



We treat this question as a problem in Bayesian model comparison flMackay||1992[ |Clark et al 



2007). The models being compared are a single gravitational wave signal present, denoted Mi, and 
no gravitational wave signals present, denoted MoQ Introduce the odds-ratio O as the the ratio of 
the probability of hypothesis Mi to the probability of the hypothesis Mo: 

O = ^^l^i (3-18a) 
p M (M |d,X) v ' 



4 Note that two or more signals present, or noise character changes, or . . ., are all different hypotheses. 
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where 

/ probability, given observations \ 

p M (M fe |d,Z) = . , ,v • tv/t • + (3-18b) 
\ a, that hypothesis is true / 

and X denotes additional, unenumerated conditions. Following Bayes' Theorem each of these prob- 
abilities can be expressed in terms of a likelihood and an appropriate a priori probability: 

PM(M!|d,X) = |r^j@ J d n & A(d\Mi,0,l)qeu(0\Mi,l) (3-19a) 

PM (M |d,X) = ^M^ A(d|M ,Z) (3-19b) 



where 



A(d|M 1) \ P r °b am lity °^ observing d assuming the gravitational \ (3 19 ) 

' ' I wave signal described by the parameters is present / 



A(d|Mo,Z) = y probability of observing d assuming no signal is present J (3-19d) 
gM(Mjt|X) = I a priori probability of hypothesis Mfc I (3-19e) 



/ a priori probability that h is described \ 

©m(0 M fc ,J) = to- i ■ i • -« t (3- 19f ) 
\ by parameters tf given hypothesis Mk I 

Z M (d\l) = ^2 PM (M k \d,l) (3-19g) 

k 

The odds-ratio O can thus be expressed as the product of two terms, one that depends only on the 
observations and one that depends only on our a priori assumptions about the outcome: 

O = B{d)p (3-20a) 

where 

B( d) = S^m^e^mu^.i) (3 20b) 

, = «<Ml£) (3.20c) 



S(d), the data-dependent contribution to 0, is referred to as the Bayes Factor (Gelman et al 



2004, pp. 184-186). The Bayes Factor reflects the evidence provided by the data d in favor of 
the hypothesis Mi relative to Mo- It is the ratio of the marginalized likelihood of the data under 
the two hypotheses Mi and Mo- When it is large compared to unity the observations favor Mi; 
when it is small compared to unity the observations favors Mo- The "odds" O are the produce of 
the Bayes factor and the priors p. Depending on our interest or prejudice p can take on different 
values. For example, if our interest is to "let the data speak for themselves" then we take p = 1; 
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i.e., we expresses no prejudice regarding the presence of absence of a gravitational wave burst in the 
data set d. Alternatively, if our interest is to express the odds in the context of some theoretical 
model or prejudice suggesting a rate of gravitational wave bursts over the period of the observation, 
then p will be a function of the rate and observation period. In any case, however, B should be 
much greater than p~ l before we are entitled to conclude with certainty that we have observed a 
gravitational wave burst. 



3.4-2. Computing the Bayes Factor 



Turn now to computing the Bayes Factor B(d) (Eq. 3-20b). Focus first on the denominator 



A(d|Mo); i.e., the probability that the particular observation d is an instance of detector network 
noise. Referring to the discussion of §3.2.1 this probability density is 



A(d|M ,X) = N{d\C) 

exp [-±d T C- l d] 



(27r) dimd det||C| 



(3-21a) 
(3-21b) 



Turn now to the Bayes Factor numerator, 



<r0A(d|Mi,0,:r) ge M(0|Mi,x), 



(3-22) 



which we recognize, upon inspection, as the evidence Zk(d\I) defined in equation (3-17b) 
The Bayes Factor is thus given by 

fZ k (d\l) 



B(d) 



A(d|M ,X) 

exp (iFC-i-dfA- 1 (RJC- 1 ^} 



d 2 n k 



i iii * ii 2dim/i_i_ 2dim/i> 

det A (j, + <7 V 



(3-23a) 
(3-23b) 



where we have taken advantage of the expression for Zh given in equation (3-13). 



3.5. Summary 

In the preceding discussion we have described a Bayesian analysis that addresses three ques- 
tions: 



1. Does the data set d include the signal from a passing gravitational wave burst? 
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2. Assuming that a gravitational wave burst is present, what is the probability that the wave is 
propagating in the direction kl 

3. Assuming a burst propagating in direction k, what is the probability that the wave at Earth 
is characterized by h? 



The answers to these questions — i.e., the principal results of this section — are given by, for the 



first question, Equation 3-23 for the second question, Equation 3-17 and, for the third question 



Equation 3-14 In the next Section we will demonstrate the effectiveness of this analysis, making 
use of these three results. 



4. Examples 
4.1. Overview 

To illustrate and demonstrate the effectiveness of the analysis techniques just described we 
apply them to simulated observations of a gravitational wave burst characteristic of the close, 
parabolic encounter of two supermassive black holes, such as might occur when the nuclear black 
holes first find each other following a major merger of two galaxies. We consider four cases: 

1. A strong signal, for which we can detect the signal, localize the source in the sky, and infer 
the radiation waveform; 

2. A moderate strength signal, for which we can detect the signal and localize the source, but 
not accurately infer the waveform; 

3. A weak signal, which can be clearly detected but not accurately localized or characterized; 
and 

4. No signal at all. 

For these examples we use the thirty pulsars in the International Pulsar Timing Array (IPTA) 
( Hobbs et al. 2009 ) as described in Table [TJ The measured timing residual for each pulsar is a 



superposition of white noise with rms timing residual given in Table [T] and red noise normalized to 
have the same spectral density as the white noise at frequency 0.2yr _1 . Of these thirty pulsars, 
ten have short-timescale timing residual noise rms less than 0.2/is, fourteen have short-timescale 
noise rms between 0.2 and 1/is and the remaining five have short-timescale noise rms between 1 
and 5/xsj^] 



5 This is a particular characterization of these pulsars based on communications at the time of this writing from the 
Parkes Pulsar Timing Array, the European Pulsar Timing Array, and the North American Nanohertz Observatory 
for Gravitational-waves. It is not a definitive characterization. We are not presenting the data associated with these 
pulsars but rather using them as an example of a realistic IPTA. 
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Table 1. The International Pulsar Timing Array pulsars, their short-timescale timing noise rms, 





and the 


telescopes 


from which those noise timing residuals were measured. 




Pulsar 


RMS (^s) 


Telescope 


1 

X 


71 Q09-3744 


054 


vl 1 J X 


2 


71 71 q_i_n747 


055 


AO 


O 


70437-471 5 
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066 


AO 


z, 

o 


71 qqq_i_9i q4 


u.uou 


vi 1 J X 


6 


7061 3-0900 


110 

U . X XVJ 


VI 1 J X 


7 


71 R4D-I-9994 


110 

U . X XVJ 


AO 


8 
o 


71 744-1 1 "34 

U X J 'X'X X X 04 


1 30 

U . X<JVJ 


ORT 

Vl 1 J X 


Q 


71 741 -Li S00 


140 


AO 


1 n 


71 600-3053 


190 


ORT 

Vl 1 J X 


1 1 

X X 




n 9nn 


AO 

JrWJ 


1 9 

X z. 


7nn' : !0-i-04 c >i 


300 

U . OVJVJ 


AO 

-Ti. vy 


1U 


TD71 1 fift^O 

JUI ll~UOOU 


^40 


X oX Ivco 


1 4 

14: 


79^1 7-Ul zl^Q 


^fiO 


AO 

JrWJ 


1 5 

X O 


1914^-07^0 


490 


X (XX iVL. O 


1 6 

XU 


71 f)1 9-1-5307 
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VJ . (Jt:VJ 


ORT 

Vl 1 J X 


1 7 


71 099-1-1 nm 

j 1UZ/ xuu x 


D 7DD 
U. t uu 


vv oxxx 


1 8 


7D91 8-1-49^9 


n 8^n 


rjRT 

Vl 1 J X 


1 Q 

X t? 


71 643-1 994 


880 

U . OOVJ 


X Cll IVL. O 


20 


J2019+2425 


0.910 


AO 


21 


J1024-0719 


0.960 


Parkes 


22 


J1455-3330 


0.960 


GBT 


23 


1918-0642 


0.960 


GBT 


24 


J1603-7202 


0.990 


Parkes 


25 
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0.990 


Parkes 


26 


J1824-2452 
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Parkes 


27 
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1.190 


Parkes 


28 


J1732-5049 


1.250 


Parkes 


29 


J1045-4509 


1.370 


Parkes 


30 


J2124-3358 


2.380 


Parkes 
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The data sets we use for these examples are constructed by 



1. Calculating the gravitational wave strain associated with the parabolic encounter of two 



supermassive black holes (see 5 4.2.1); 



2. Evaluating the gravitational wave contribution to the pulse arrival time for each pulsar de- 
scribed in Table [TJ 



3. Adding the appropriate noise to the "gravitational wave" timing residuals (see 54.2.2); 

4. Removing the best-fit linear trend from the noisy timing residuals. 



At present, actual pulsar timing residual observations are constructed by fitting actual pulse time- 
of-arrival data for each pulsar to a timing model characterized by, among other parameters, the 



pulsar period and period derivative (Edwards et al. 2006). The final step in the construction of 
our simulated data — removing the linear trend — modifies the data in a manner similar to the 
"fitting-out" procedure that occurs in the construction of actual timing residual data sets. 

To summarize, our simulated data sets model — in schematic form — the major features of 
modern pulsar timing array data sets and the elements that complicate their analysis: white timing 
noise on short timescales, red timing noise on long timescales, and formation of timing residuals 
through fitting pulse arrival times to a global timing model. 



4.2. Construction of simulated data sets 

4-2.1. Parabolic encounter of two supermassive black holes 

Following the major merger of two galaxies, each harboring a nuclear supermassive black hole, 
dynamical friction will drive the nuclear black holes to the nucleus of the merged galaxy. Eventually 
they will find each other, form a binary, and coalesce. When they first find each other there may 
occur a series of close, high-speed encounters, each leading to a burst of radiation, whose duration 
may be estimated as twice the ratio of the impact parameter to the velocity at periapsis. We 
adopt this burst as an exemplar for the purpose of demonstrating the effectiveness of the analysis 
techniques just described gravitational wave burst. At the same time, however, we emphasize 
that the parabolic encounter gravitational wave model used here is intended as a stand-in for any 
gravitational wave burst: i.e., the particular model and model parameters adopted here do not 
correspond to a case we regard as realistic. 

We model the parabolic encounter radiation burst via the quadrupole formula applied to the 
Keplerian parabolic trajectories of the equivalent Newtonian system. In the quadrupole approxi- 
mation the gravitational waves radiated near periapsis are projections of the second time derivative 
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of the systems quadrupole moment: i.e., 



h+ 



(4- la) 
(4-lb) 



where k is the unit vector in the direction of wave propagation, we have adopted the Einstein 
summation convention of summing over repeated indices, and work in units where G = c = 1. For 
Keplerian parabolic orbits the trajectories (and, correspondingly, the system's quadrupole moment) 
can be expressed in closed form. Without loss of generality we take the system's motion to be in 
the xy plane and the periapsis at y = and x > 0, in which case 

qxx = \-3w\ ( w f - 6u;? + 24wl - 16) + w (7wf - 30wf + 24wf - 16)1 (4-2a) 

WqW{0 



qvv 
Q 



[Swf (wj -A)+w (5wi - 4)] 



WQivfb 



" 3 -i6>/2 



[w (-I8wf + 32wj + 32) + 3wi (7wf - 30wf + 2Awf + 16)] 



(4,2b) 
(4-2c) 



WqWiI 

where M and \i are the system's total and reduced mass, b is the impact parameters, and ui and 
wi are given by the following: 



WQ 



+ 9- 



M ft 
~b 



t M 

3 bVT + Wo 



1/3 



(4-2d) 
(4-2e) 



Similarly, the gravitational wave contribution to the timing residual is a projection of the time inte- 
gral of h (see Eq. 2-5), which is proportional to the first time derivative of the system's quadrupole 
moment: 



QXX 

qvv 
Qxy 



/ib 



\^2wowf V b 
T7 

(«•,' I) 

'M 



(wf - 4) (wf - 6wj + 4) 



wqw\ V b 



V2wowf 




3wf + 8wf + 16wj - 24) 



(4-3a) 
(4-3b) 
(4-3c) 



4-2.2. Timing noise 



The millisecond pulsars used in modern pulsar timing arrays typically show white timing noise 
on short timescales, turning to red noise on timescales of 5-10 years. For the demonstrations 
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here we model the timing noise as the superposition of white noise and red noise, with the red 
noise contribution normalized to have the same amplitude as the white noise contribution at the 
frequency / re< j = 0.2 yr -1 . With this normalization the noise power spectrum for each pulsar in 
our array is completely determined by the short-timescale (white) timing noise rms given in Table 

m 

To compute the red contribution to the timing noise we applying a digital integrator to white 



noise. To design the integrator we follow Tseng (2006), choosing a single sub-division of the unit 



delay, a seventh-order FIR filter, and a cascade of three unit delays. The corresponding integrator 
is given by the transfer function 

1 -5 + 49z _1 - 245z~ 2 + 1225z~ 3 + 1225z" 4 - 245z" 5 + 49z" 6 - 5z~ 7 , A A . 

H{Z) = 2P 6(1^) (4 " 4) 

Figure [l] shows the characteristics of the power spectral density of the simulated timing noise 
normalized for PSR J1909-3744. 



4.3. Analysis of simulated data for strong, moderate, weak, and no signal 

Our analysis methodology is designed to answer three questions: (1) is a gravitational wave 
signal present? (2) where is the source? and (3) what is the detailed structure of the waveform? 
Here we explore the ability of our analysis to answer these questions. For weak signals it may be 
possible to answer definitively the first of these questions, while being unable to answer the second 
or third. For stronger signals it may be possible to answer the first question definitively, the second 
moderately well, and the third not at all. Finally, for the strongest signals all three questions may 
be answered in detail. We illustrate all three cases in the three following subsections, beginning 
with a strong signal example and ending with a weak signal example. In each case our "source" 
has waveform characteristic of the parabolic encounter of two 10 9 M Q black holes, with impact 
parameter 180 M (0.02 pc), orbital plane face-on to the Earth line-of-sight, and in the direction of 
the Virgo cluster (RA 12. 5h, dec 12.5deg). Figure [2] shows, in two panels, the gravitational wave 
strain incident at Earth (top panel) and the induced timing residuals in a sample of six of the thirty 
IPTA pulsars when the source is at a distance of 15 Mpc. We conclude with a subsection exploring 
how the analysis performs when applied to a data set containing no signal at all. 



4-3.1. Strong signal 

Figures [2] through [5] and the first row of Table [2] summarize the results of applying the method- 
ology described in Section || to a pulsar timing array dataset including a strong "flyby" signal, 
constructed as described in Section §4.2| In this strong signal case the source is placed at a dis- 
tance of 15 Mpc. The top panel of Figure [2] shows the strain incident at Earth and the bottom 
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Fig. 1. — Power spectral density of simulated timing noise for PSR J1909-3744. The noise is 
simulated as the sum of a white noise contribution, which is determining at high frequencies, and a 
red noise contribution, which is determining at low frequencies. The cross-over frequency is chosen 
to be 0.2 yr _1 , which is characteristic of timing array millisecond pulsars. For more details see 
pX2 
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panel the timing residual induced in a sample of six of the thirty IPTA pulsars. Figure [3] shows the 
same timing residuals, from the same selection of pulsars, embedded in the "red-plus-white" timing 



noise described in {4.2.2 For this strong signal the gravitational wave induced timing residuals are 
readily apparent in the quietest of the IPTA pulsars (e.g., top two panels of Fig. [5]), less so in the 
pulsars with moderate timing noise (middle panels of Fig. [3]), and much less so in the pulsars with 
large timing noise (the bottom two panels of Fig. [3]). 



Applying the analysis described in Section 3.4 to this "strong signal" dataset instance we find 



(for our particular instantiation of noise) that the Bayes Factor has a value of exp(3.8 x 10 3 ), 
corresponding to overwhelming evidence for the presence of a gravitational wave in this data set. 



Having concluded that a signal is present we use the analysis described in Section 3.3 to localize 
the source. Figure [4] shows the results of this analysis as the natural log of the probability density 
that the source is in the direction f2. Also shown is the smallest contour containing 90% of the 
total probability, whose area is much less than 1 deg 2 . 

Finally, having detected the source and localized it on the sky, we apply the analysis of Section 
3.2 to infer the radiation waveform. Figure [5] shows the result of this analysis, superposed with the 



actual radiation waveform. In this example the gravitational wave strain is identified with a power 
signal-to- noise of 38. 



4-3.2. Moderate signal 

Figures [6] through [8] show the results of applying the methodology described in Section [3] to a 
moderate strength "flyby" signal observed in the current IPTA. In this case the source is placed at 
a distance of 100 Mpc in the direction of the Virgo Cluster. Figure [2j with the appropriate scaling 
of the abscissae (i.e., by 15 mpc/100 mpc) shows the gravitational wave strain incident on the IPTA 
and the corresponding induced timing residuals in a selection of IPTA pulsars. Figure [6] shows the 
timing residuals, from the same selection of pulsars as in Figure [2| embedded in the red-plus- white 



timing noise described in {4.2.2 For this moderate strength signal the gravitational wave induced 
timing residuals are apparent in the quietest of the IPTA pulsars (e.g., top two panels of Fig. [6|, 
but not apparent in the residuals with of the other pulsars. 



Applying the analysis described in Section 3.4 to this data set we find that the Bayes Factor has 



a value of exp(66.), again corresponding to overwhelming evidence for the presence of a gravitational 
wave signal. 

Having concluded that a signal is present we attempt to localize the source using the analysis 
described in Section 3.3 Figure [7] shows the results of our localization analysis as the log of the 



probability density that the source is in the direction Q. Contours enclosing the smallest area 
containing 90% of the total probability are also shown. These contours, which encloses an area of 
5.8 x 10 2 deg 2 , correctly include the actual source location. 
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Fig. 2. — The gravitational wave strain incident on Earth and the corresponding timing residuals 
induced in a sample of International Pulsar Timing Array pulsars. In this example the waves are 
characteristic of the parabolic encounter of two 10 9 M Q black holes, impact parameter 180 M (i.e., 
0.02 pc), at a distance of 15 Mpc in the direction of the Virgo Cluster of galaxies (RA 12h5m, 



Dec 12.5deg). (See discussion of £4.3.1 ) The top panel shows the radiation waveform in the two 
independent polarization states. The bottom panel shows the timing residuals induced by the 
waveform in a sample of 6 of the 30 IPTA pulsars. 




Fig. 3. — The superposition of the gravitational wave induced timing residuals, for the same sample 
of IPTA pulsars shown in the bottom panel of Fig. [2j with the "red+white" timing noise (see { 4.2.2 ) 
characteristic of typical millisecond pulsar timing noise. 
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Fig. 4. — Natural log of the inferred probability density that the source of gravitational waves 
present in the "strong signal" simulated IPTA data set described in £4.3.1 is found at location f2. 
The smallest 90% probability contour has an area much less than 1 deg 2 and includes the actual 
source location. The white squares show the locations of the thirty IPTA pulsars used as detectors. 
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Table 2. Results summary for data analysis applied to four simulated datasets. In all cases the 
signal corresponds to radiation from a parabolic fly-by of two 10 9 Mq black holes propagating 
from the direction of the Virgo Cluster. In the "strong" signal case the source is located at 
15 Mpc; in the "moderate" signal case the source is at a distance of 100 Mpc; in the "weak" 
signal case the source is at a distance of 260 Mpc; and in the final, "absent" signal case the 



simulated data set consists of timing noise alone. For details see section {4.3 



Signal In 5(d) Af2 90 % (deg 2 ) p 2 a+ a x 



Strong 


3.8 


X 


10 3 




<C 1 


3.8 


x 


io- 


-i 


1.1 


X 


io- 


-13 


1.2 


X 


io- 


-13 


Moderate 


6.6 


X 


10 1 


5.8 


x 10 2 


8.7 


X 


10" 


-i 


2.4 


X 


io- 


-14 


2.4 


X 


io- 


-14 


Weak 


2.2 


X 


10° 


4.2 


x 10 3 


2.1 


X 


10" 


-i 


2.0 


X 


io- 


-14 


2.0 


X 


io- 


-14 


Absent 


-8.4 


X 


10° 


1.2 


x 10 4 


7.9 


X 


io- 


-2 


1.8 


X 


io- 


-14 


1.9 


X 


io- 


-11 





Fig. 5. — The inferred h+ and h x radiation waveforms for the example data set described in {■ 



4.3.1 



- 28 - 



Finally, having detected the source and localized it on the sky, we apply the analysis of Section 
3.2 to infer the radiation waveform. Figure [8] shows the result of this analysis, made assuming we 
know the actual source location, superposed with the actual radiation waveform. In this example 
the power signal-to-noise ratio is 0.87, which confirms the impression given by the figure that this 
inference is not significant. 

This moderate signal amplitude case shows clearly a regime where the gravitational wave burst 
is strong enough to be unambiguously detected and the general direction to the source clearly 
identified (even if not so precisely that an optical counterpart may be sought), but not strong 
enough to characterize the burst waveform. 



4-3.3. Weak signal 

Finally, we consider a dataset that includes a signal at the edge of detectability: i.e., a data 
set where the Bayes Factor corresponds to 9:1 odds of a signal being present, with results shown in 



Figures |9] through 11 In this case our binary source is placed at a distance of 2.6 x 10 2 Mpc. Figure 
[2J with the appropriate scaling of the abscissae (i.e., by 15 kpc/261 mpc) shows the gravitational 
wave strain incident on the IPTA and the corresponding induced timing residuals in a selection of 
IPTA pulsars. Figure [9] shows the timing residuals from the same selection of pulsars as in Figure 
[2j embedded in white timing noise with rms given in Table [T] For this weak signal the gravitational 
wave induced timing residuals are not readily apparent even in the quietest pulsars (e.g., top panel 
of Fig.|9b. 



Applying the analysis described in Section 3.4 to this data set we find that the Bayes Factor has 



a value of exp(2.2) = 9. At this level our prejudice regarding the likelihood of gravitational wave 
bursts passing through our timing array plays a critical role in deciding whether the overall odds 
- i..e, the product of the Bayes Factor with the "expectation odds" — are in favor of detection or 
not. Supposing that they are we next attempt to localize the source using the analysis described 



in Section 3.3 Figure 10 shows the results of our localization analysis as the log of the probability 
density that the source is in the direction O. In this case, the 90% contour encloses an area of 
4.2 x 10 3 deg 2 scattered about the sky: i.e., the wave, while strong enough to be detected, is not 
strong enough to be localized. 



Finally, and for completeness, we apply the analysis of Section 3.2 to infer the radiation 



waveform. Figure 11 shows the result of this analysis, made assuming that we know the actual 
source location on the sky, superposed with the actual radiation waveform. In this example the 
power signal-to- noise ratio is 0.21. 
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Fig. 6. — The superposition of the gravitational wave induced timing residuals, for the same sample 
of IPTA pulsars shown in the bottom panel of Fig. [2j with the "red+white" timing noise (see { 4.2.2 ) 
characteristic of typical millisecond pulsar timing noise. 
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Fig. 7. — Natural log of the inferred probability density that the source of gravitational waves 
present in the "moderate signal" simulated IPTA data set described in £4.3.2 is found at location 
0,. Also shown is the smallest 90% probability contour, which has an area of 5.8 x 10 2 deg 2 . The 
white squares show the locations of the thirty IPTA pulsar baselines used as detectors. 
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Fig. 8. — The inferred h+ and h x radiation waveforms for the example data set described in £4.3.2 
While strong enough to be detected the signal is too weak for us to infer its waveform. 
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Fig. 9. — The superposition of the gravitational wave induced timing residuals, for the same sample 
of IPTA pulsars shown in the bottom panel of Fig. [2j with the "red+white" timing noise (see { 4.2.2 ) 
characteristic of typical millisecond pulsar timing noise. 
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Fig. 10. — Natural log of the inferred probability density that the source of gravitational waves 
present in the "weak signal" simulated IPTA data set described in §4.3.3 is found at location 
f2. While strong enough to be detected, the signal is too weak to be reliably localized: the 90% 
probability contour has an area of 4.2 x 3 deg 2 . The white squares show the locations of the thirty 
IPTA pulsar baselines used as detectors. 
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Fig. 11. — The inferred h + and h x radiation waveforms for the example data set described in 
£4.3.3 While strong enough to be detected, the signal is too weak to for us to infer its waveform. 
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4-3.4- No signal 

Finally, we apply our analysis to a data set consisting of noise only. In this particular instance 
the Bayes Factor is exp(— 8.4) = 2.2 x 10 -4 : i.e., overwhelming evidence for the absence of a 
gravitational wave burst. For completeness, under the assumption that a single source is present 



we show in Figure 12 the inferred probability density for the source location on the sky. In this case 
the 90% confidence interval has an area of 1.2 x 10 4 deg 2 : i.e., approximately 1/3 of the sky. Lastly, 
under the assumption that there is a source in the direction of the Virgo Cluster we attempt to 



infer a waveform from this data set. Figure 13 shows the results of this analysis, which correspond 
to a power signal-to-noise ratio of 7.8 x 1CP 2 . We conclude that the analysis described here is fully 
capable of identifying data sets that contain no evidence of a gravitational wave signal. 



5. Conclusions 



In the history of astronomy few (if any) new observational windows have been as eagerly 
anticipated as the gravitational wave window, whose opening will provide us with a novel and direct 
view of astronomical phenomena that can now be inferred at best dimly and indirectly. Making 
sense of what we see through this new window requires analysis tools and techniques adapted to 
the unique nature of our new "telescopes" and the sources they enable us to study. Here we have 
described an intra-related suite of analysis techniques for gravitational wave astronomy designed 
to address quantitatively three specific questions: 

1. What are the "odds" that a gravitational wave detector data set includes the signal from a 
gravitational wave burst? 

2. Assuming that a gravitational wave burst is present in a data set, what is the probability 
that the wave is propagating in direction fe? 

3. Assuming the presence of a burst propagating in direction k, what is the probability that the 
wave at Earth is characterized by the functions h+(u) and h x (u), u = t — k ■ x, representing 
the + and x polarization state waveforms? 



We address these questions in the specific context of gravitational wave detection using pulsar 
timing array data. Until recently, analyses for gravitational wave detection using timing data from 
an array of pulsars has focused on stationary sources: e.g., a stochastic gravitational wave signal 
or the signal from a binary system. By addressing burst sources we also add to the very recent 



literature examining how pulsar timing data can be used to detect gravitational wave bursts (van 



Haasteren & Levin 2009) such as might arise a close fly-by or collision of two supermassive black 



holes or from a cosmic string cusp (Binetruy et al. 2009 Key & Cornish 2009). 



To demonstrate the efficacy of our analysis we applied it to four synthetic timing residual data 



sets representative of observations using the International Pulsar Timing Array (Hobbs et al. 112009 



- 36 - 




Fig. 12. — Natural log of the inferred probability density that a source of gravitational waves is 
present in the direction f2 for a data set consisting only of noise. The 90% probability contour has 
an area of 1.2 x 4 deg 2 . The white squares show the locations of the thirty IPTA pulsar baselines 
used as detectors. See 94.3.41 for more details. 
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Fig. 13. — The inferred h+ and h x radiation waveforms for a data set consisting only of timing 
noise. The inferred waveform corresponds to a signal-to-noise of 7.8 x 1(T 2 See UIm for further 
discussion. 
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IPTA). Each data set included simulated timing noise, constructed to be characteristic of actual 
IPTA timing noise. Three of the datasets included the timing residual signature of a gravitational 
wave burst characteristic of a parabolic fly-by of two supermassive black holes; the fourth did not. 
The three "signal-present" cases varied only by the gravitational wave signal amplitude. In the case 
of the strongest signal the burst was unambiguously detected, localized to much better than a deg 2 , 
and the waveform in the individual polarization states recovered. In the moderate signal amplitude 
case the signal was, again, unambiguously detected and the general direction to the source clearly 
determined; however, the signal amplitude was too low to infer the waveform characteristics. In 
the third case the signal was strong enough to be detected but too weak to be characterized or to 
allow the source to be localized. Finally, in analyzing noise alone the calculated odds were, as they 
should be, unambiguously against the presence of a gravitational wave burst. 

At present, pulsar timing array data sets are constructed by fitting a timing mode to the time of 
arrival (TOA) data for each pulsar. This timing model includes, in parameterized form, all the non- 
gravitational- wave contributions that affect the pulse arrival times. The residual differences between 
the timing model predictions and the actual arrival times are then analyzed for the signature of a 
passing gravitational wave. This procedure has the disadvantage of being incapable of identifying 
any gravitational wave whose effect on the arrival time of individual pulsars is degenerate with 
any of the non-gravitational-wave effects that are part of the timing model. Our analysis may 
be extended to infer, in addition to the gravitational radiation waveform, the other timing model 
parameters. This extended analysis may be applied to TOA data directly, avoiding entirely the 
problem of "fitting-out" gravitational wave contributions whose character is similar to other timing 
model contributions. We intend to investigate this extension in future work. 

While our presentation and discussion have focused on pulsar timing array observations the 
analysis methodology that we describe applies equally well and without modification to gravitational 



wave data taken from ground-based detector networks like the LIGO- Virgo network (Accadia et al. 



20101 |Riles et al.||2010| ), or for the analysis of LISA QJennrich||2009t |Merkowitz et~al^|2009| data. 



We thank Lynn Baker for valuable conversations, and Martin Hendry and Graham Woan for 
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Foundation grants AST-0748580 (AL) and PHY 06-53462 (LSF). 



REFERENCES 

Accadia, T., Acernese, F., Antonucci, F., Aoudia, S., Arun, K. G., Astone, P., Ballardin, G., 
Barone, F., Barsuglia, M., Bauer, T. S., Beker, M. G., Belletoile, A., Bigotta, S., Birindelli, 
S., Bizouard, M. A., Blom, M., Boccara, C, Bondu, F., Bonelli, L., Bonnand, R., Bosi, L., 
Braccini, S., Bradaschia, C, Brillet, A., Brisson, V., Budzyhski, R., Bulik, T., Bulten, H. J., 
Buskulic, D., Buy, C, Cagnoli, G., Calloni, E., Campagna, E., Canuel, B., Carbognani, 
F., Cavalier, F., Cavalieri, R., Cella, G., Cesarini, E., Chassande-Mottin, E., Chincarini, 



-39- 



A., Cleva, F., Coccia, E., Colacino, C. N., Colas, J., Colla, A., Colombini, M., Corsi, A., 
Coulon, J., Cuoco, E., D'Antonio, S., Dari, A., Dattilo, V., Davier, M., Day, R., De Rosa, 
R., del Prete, M., Di Fiore, L., Di Lieto, A., Emilio, M. D. P., Di Virgilio, A., Dietz, A., 
Drago, M., Fafone, V., Ferrante, I., Fidecaro, F., Fiori, I., Flaminio, R., Fournier, J., Franc, 
J., Frasca, S., Frasconi, F., Freise, A., Galimberti, M., Gammaitoni, L., Garufi, F., Gemme, 
G., Genin, E., Gennai, A., Giazotto, A., Gouaty, R., Granata, M., Greverie, C., Guidi, G., 
Heitmann, H., Hello, P., Hild, S., Huet, D., Jaranowski, P., Kowalska, I., Krolak, A., La 
Penna, P., Leroy, N., Letendre, N., Li, T. G. F., Lorenzini, M., Loriette, V., Losurdo, G., 
Mackowski, J., Majorana, E., Man, N., Mantovani, M., Marchesoni, F., Marion, F., Marque, 
J., Martelli, F., Masserot, A., Menzinger, F., Michel, C, Milano, L., Minenkov, Y., Mohan, 
M., Moreau, J., Morgado, N., Morgia, A., Mosca, S., Moscatelli, V., Mours, B., Neri, I., 
Nocera, F., Pagliaroli, G., Palomba, C, Paoletti, F., Pardi, S., Parisi, M., Pasqualetti, 
A., Passaquieti, R., Passuello, D., Persichetti, G., Pichot, M., Piergiovanni, F., Pietka, M., 
Pinard, L., Poggiani, R., Prato, M., Prodi, G. A., Punturo, M., Puppo, P., Rabaste, O., 
Rabeling, D. S., Rapagnani, P., Re, V., Regimbau, T., Ricci, F., Robinet, F., Rocchi, A., 
Rolland, L., Romano, R., Rosihska, D., Ruggi, P., Sassolas, B., Sentenac, D., Sturani, R., 
Swinkels, B., Toncelli, A., Tonelli, M., Tournefier, E., Travasso, F., Trummer, J., Vajente, 
G., van den Brand, J. F. J., van der Putten, S., Vavoulidis, M., Vedovato, G., Verkindt, D., 
Vetrano, F., Vicere, A., Vinet, J., Vocca, H., Was, M., & Yvert, M. 2010, Journal of Physics 
Conference Series, 203, 012074 

Anholm, M., Ballmer, S., Creighton, J. D. E., Price, L. R., & Siemens, X. 2009, Phys. Rev. D, 79, 
084030 

Astone, P., Baggio, L., Bassan, M., Bignotto, M., Bonaldi, M., Bonifazi, P., Cerdonio, M., Coccia, 
E., Conti, L., D'Antonio, S., Emilio, M. d. P., Drago, M., Fafone, V., Falferi, P., Foffa, S., 
Fortini, P., Frasca, S., Giordano, G., Hamilton, W. O., Hanson, J., Johnson, W. W., Liguori, 
N., Longo, S., Maggiore, M., Marin, F., Marini, A., McHugh, M. P., Mezzena, R., Miller, 
P., Minenkov, Y., Mion, A., Modestino, G., Moleti, A., Nettles, D., Ortolan, A., Pallottino, 
G. V., Pizzella, G., Poggi, S., Prodi, G. A., Re, V., Rocchi, A., Ronga, F., Salemi, F., 
Sturani, R., Taffarello, L., Terenzi, R., Vedovato, G., Vinante, A., Visco, M., Vitale, S., 
Weaver, J., Zendri, J. P., & Zhang, P. 2010, ArXiv e-prints 

Binetruy, P., Bohe, A., Hertog, T., & Steer, D. A. 2009, Phys. Rev. D, 80, 123510 

Bondarescu, R., kumar Kopparpu, R., & Finn, L. S. 2010, Estimating gravitational wave Stokes 
Parameters with ground-based detectors, in preparation 

Cawley, G. C. k Talbot, N. L. C. 2007, Journal of Machine Learning Research, 8, 841 

Christodoulou, D. 1991, Phys. Rev. Lett., 67, 1486 

Clark, J., Heng, I. S., Pitkin, M., & Woan, G. 2007, Phys. Rev. D, 76, 043003 



-40- 



Damour, T. & Vilenkin, A. 2001, Phys. Rev. D, 64, 064008 
Demorest, P. B. 2007, PhD thesis, University of California, Berkeley 
Detweiler, S. 1979, ApJ, 234, 1100 

Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549 

Finn, L. S. 2009, Phys. Rev. D, 79 

Foster, R. S. & Backer, D. C. 1990, ApJ, 361, 300 

Galatsanos, N. P. & Katsaggelos, A. K. 1992, IEEE Transactions on Image Processing, 1, 322 

Galatsanos, N. P., Mesarocic, V. Z., Molina, R., & Katsaggelos, A. K. 1998, in Proc. SPIE, Vol. 
3459, Bayesian Inference for Inverse Problems, ed. A. M. Djafari, 337-348 

Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. 2004, Bayesian Data Analysis, 2nd edn., 
Texts in Statistical Science (Boca Raton, Florida: Chapman & Hall/CRC) 

Hellings, R. W. & Downs, G. S. 1983, ApJL, 265, L39 

Hobbs, G., Archibald, A., Arzoumanian, Z., Backer, D., Bailes, M., Bhat, N. D. R., Burgay, 
M., Burke-Spolaor, S., Champion, D., Cognard, I., Coles, W., Cordes, J., Demorest, P., 
Desvignes, G., Ferdman, R. D., Finn, L., Freire, P., Gonzalez, M., Hessels, J., Hotan, A., 
Janssen, G., Jenet, F., Jessner, A., Jordan, C, Kaspi, V., Kramer, M., Kondratiev, V., 
Lazio, J., Lazaridis, K., Lee, K. J., Levin, Y., Lommen, A., Lorimer, D., Lynch, R., Lyne, 
A., Manchester, R., McLaughlin, M., Nice, D., Oslowski, S., Pilia, M., Possenti, A., Purver, 
M., Ransom, S., Reynolds, J., Sanidas, S., Sarkissian, J., Sesana, A., Shannon, R., Siemens, 
X., Stairs, I., Stappers, B., Stinebring, D., Theureau, G., van Haasteren, R., van Straten, 
W., Verbiest, J. P. W., Yardley, D. R. B., & You, X. P. 2009, ArXiv e-prints 

Hobbs, G., Jenet, F., Lommen, A., Coles, W., Verbiest, J. P. W., k Manchester, R. 2008, in 
American Institute of Physics Conference Series, Vol. 983, 40 Years of Pulsars: Millisecond 
Pulsars, Magnetars and More, ed. C. Bassa, Z. Wang, A. Cumming, & V. M. Kaspi, 630-632 

Jaffe, A. H. & Backer, D. C. 2003, ApJ, 583, 616 

Jenet, F., Hobbs, G., Lee, K., & Manchester, R. 2005, ApJ, 625, L123 
Jenet, F. A., Creighton, T., & Lommen, A. 2005a, ApJL, 627, L125 

Jenet, F. A., Hobbs, G. B., van Straten, W., Manchester, R. N., Bailes, M., Verbiest, J. P. W., 
Edwards, R. T., Hotan, A. W., Sarkissian, J. M., & Ord, S. M. 2006, ApJ, 653, 1571 

Jenet, F. A., Lommen, A., Larson, S. L., & Wen, L. 2004, ApJ, 606, 799 



- 41 - 



Jenet, F. A., Lommen, A., Larson, S. L., & Wen, L. 2005b, in Astronomical Society of the Pacific 
Conference Series, Vol. 328, Binary Radio Pulsars, ed. F. A. Rasio & I. H. Stairs, 399 — h 

Jennrich, O. 2009, Class. Quantum Grav., 26, 153001 

Kawamura, S., Ando, M., Nakamura, T., Tsubono, K., Tanaka, T., Funaki, I., Seto, N., Numata, 
K., Sato, S., Ioka, K., Kanda, N., Takashima, T., Agatsuma, K., Akutsu, T., Akutsu, T., 
Aoyanagi, K., Arai, K., Arase, Y., Araya, A., Asada, H., Aso, Y., Chiba, T., Ebisuzaki, T., 
Enoki, M., Eriguchi, Y., Fujimoto, M., Fujita, R., Fukushima, M., Futamase, T., Ganzu, 
K., Harada, T., Hashimoto, T., Hayama, K., Hikida, W., Himemoto, Y., Hirabayashi, H., 
Hiramatsu, T., Hong, F., Horisawa, H., Hosokawa, M., Ichiki, K., Ikegami, T., Inoue, K. T., 
Ishidoshiro, K., Ishihara, H., Ishikawa, T., Ishizaki, H., Ito, H., Itoh, Y., Kamagasako, 
S., Kawashima, N., Kawazoe, F., Kirihara, H., Kishimoto, N., Kiuchi, K., Kobayashi, S., 
Kohri, K., Koizumi, H., Kojima, Y., Kokeyama, K., Kokuyama, W., Kotake, K., Kozai, 
Y., Kudoh, H., Kunimori, H., Kuninaka, H., Kuroda, K., Maeda, K., Matsuhara, H., Mino, 
Y., Miyakawa, O., Miyoki, S., Morimoto, M. Y., Morioka, T., Morisawa, T., Moriwaki, S., 
Mukohyama, S., Musha, M., Nagano, S., Naito, I., Nakagawa, N., Nakamura, K., Nakano, 
H., Nakao, K., Nakasuka, S., Nakayama, Y., Nishida, E., Nishiyama, K., Nishizawa, A., 
Niwa, Y., Ohashi, M., Ohishi, N., Ohkawa, M., Okutomi, A., Onozato, K., Oohara, K., 
Sago, N., Saijo, M., Sakagami, M., Sakai, S., Sakata, S., Sasaki, M., Sato, T., Shibata, M., 
Shinkai, H., Somiya, K., Sotani, H., Sugiyama, N., Suwa, Y., Tagoshi, H., Takahashi, K., 
Takahashi, K., Takahashi, T., Takahashi, H., Takahashi, R., Takahashi, R., Takamori, A., 
Takano, T., Taniguchi, K., Taruya, A., Tashiro, H., Tokuda, M., Tokunari, M., Toyoshima, 
M., Tsujikawa, S., Tsunesada, Y., Ueda, K., Utashima, M., Yamakawa, H., Yamamoto, K., 
Yamazaki, T., Yokoyama, J., Yoo, C, Yoshida, S., k, Yoshino, T. 2008, Journal of Physics 
Conference Series, 122, 012006 

Keren, D. & Werman, M. 1996, in Maximum Entropy and Bayesian Methods, ed. K. M. Hanson & 
R. N. Silver (Netherlands: Kluwer Academic Publishers), 77-84 

Key, J. S. & Cornish, N. J. 2009, Phys. Rev. D, 79, 043014 

Kittel, C. 1958, Elementary Statistical Physics (New York: Wiley) 

Kuroda, K. k the LCGT Collaboration. 2006, Class. Quantum Grav., 23, 215 

Leblond, L., Shlaer, B., & Siemens, X. 2009, Phys. Rev. D, 79, 123519 

Lommen, A. N. 2001, PhD thesis, University of California, Berkeley, Berkeley, CA, U.S.A. 

Lommen, A. N. & Backer, D. C. 2001, ApJ, 562, 297 

Lommen, A. N., Backer, D. C, Splaver, E. M., & Nice, D. J. 2003, in Astronomical Society of 
the Pacific Conference Series, Vol. 302, Radio Pulsars, ed. M. Bailes, D. J. Nice, & S. E. 
Thorsett, 81-+ 



- 42 - 



Mackay, D. J. C. 1992, Neural Computation, 4, 415 

MacKay, D. J. C. 1996, in Maximum Entropy and Bayesian Methods: Santa Barbara, California, 
USA, 1993, ed. G. R. Heidbreder, Fundamental Theories of Physics (Springer), 41 

MacKay, D. J. C. 1999, Neural Computation, 11, 1035 

McHugh, M. P., Zalamansky, C, Vernotte, F., k Lantz, E. 1996, Phys. Rev. D, 54, 5993 

Merkowitz, S. M., Castellucci, K. E., Depalo, S. V., Generie, J. A., Maghami, P. G., k Peabody, 
H. L. 2009, Journal of Physics Conference Series, 154, 012021 

Misner, C. W., Thorne, K. S., k Wheeler, J. A. 1973, Gravitation (San Francisco: W. H. Freeman 
k Co.) 

Molina, R., Katsaggelos, A. K., k Mateos, J. 1999, IEEE Transactions on Image Processing, 8, 231 

Riles, K., the LIGO Scientific Collaboration, k the Virgo Scientific Collaborations. 2010, Journal 
of Physics Conference Series, 203, 012002 

Saulson, P. 1994, Fundamentals of Interferometric Gravitational Wave Detectors (Signapore: World 
Scientific) 

Sazhin, M. V. 1978, Soviet Astronomy, 22, 36 

Sesana, A., Vecchio, A., k Colacino, C. N. 2008, MNRAS, 390, 192 
Siemens, X., Mandic, V., k Creighton, J. 2007, Phys. Rev. Lett., 98, 111101 

Smith, J. R. k the LIGO Scientific Collaboration. 2009, Class. Quantum Grav., 26, Numerical 
Relativity Data Analysis Meeting, Syracuse, NY, MAY 11-AUG 14, 2008 

Stebbins, R. 2006, in AIP Conference Proceedings, Vol. 873, Laser Interferometer Space Antenna: 
6th International LISA Symposium, ed. S. M. Merkowitz k J. C. Livas (American Institute 
of Physics), 3-12 

Summerscales, T. Z., Burrows, A., Finn, L. S., k Ott, C. D. 2008, ApJ, 678, 1142 

Thompson, A. M. k Kay, J. 1993, Inverse Problems, 9, 749 

Thorne, K. S. 1992, Phys. Rev. D, 45, 520 

Thorsett, S. E. k Dewey, R. J. 1996, Phys. Rev. D, 53, 3468 

Tseng, C.-C. 2006, IEE Proceedings- Vision, Image and Signal Processing, 153, 79 

van Haasteren, R. k Levin, Y. 2009, ArXiv e-prints 

van Haasteren, R., Levin, Y., McDonald, P., k Lu, T. 2009, MNRAS, 395, 1005 



Wiseman, A. G. & Will, C. M. 1991, Phys. Rev. D, 44, R2945 
Wyithe, J. S. B. & Loeb, A. 2003, ApJ, 590, 691 



This preprint was prepared with the AAS M^X macros v5.2. 



